function [] = runCombineCounterfactuals(runId, cnfactTypes, plotName, varargin)

if length(varargin) >= 1
	includeBase = varargin{1};
else
	includeBase = true;
end

%%% Add paths
addpath('../general_funcs');
addpath('../consumerAdoptions');
addpath('../providerAdoptions');
addpath('../providerExits');

%%% This script loads the results of dynamic prediction in the base case and in a set of counterfactual scenarios.
%%% It then produces time series plots to investigate the effect of interventions / what would have happened in alternative scenarios.

%%% Get input / output paths
config = getConfig_jointProcess(runId);

basecaseFolder = sprintf('%s/dynamicPredictions', config.runFolder);

cnfactNames = sprintfc('cnfact%d', cnfactTypes);

outputFolder = createFolderIfNotExist(sprintf('%s/cnfacts/plots/%s', config.runFolder, plotName));

NumCnfacts = length(cnfactNames);

%%% Load base-case results
load(sprintf('%s/pred_vals.mat', basecaseFolder), 'pred_aggVals', 'periodLabels');

M.consumerAdoptions = cell2table(periodLabels);
M.providerAdoptions  = cell2table(periodLabels);
M.providerExits      = cell2table(periodLabels);
M.laggedConsumers   = cell2table(periodLabels);
M.laggedProviders    = cell2table(periodLabels);

if includeBase
	M2.consumerAdoptions = array2table(mean(pred_aggVals.consumerAdoptions, 2)); M2.consumerAdoptions.Properties.VariableNames = {'base'};
	M2.providerAdoptions  = array2table(mean(pred_aggVals.providerAdoptions,  2)); M2.providerAdoptions.Properties.VariableNames  = {'base'};
	M2.providerExits      = array2table(mean(pred_aggVals.providerExits,      2)); M2.providerExits.Properties.VariableNames      = {'base'};
	M2.laggedConsumers   = array2table(mean(pred_aggVals.laggedConsumers,   2)); M2.laggedConsumers.Properties.VariableNames   = {'base'};
	M2.laggedProviders    = array2table(mean(pred_aggVals.laggedProviders,    2)); M2.laggedProviders.Properties.VariableNames    = {'base'};
	
	M.consumerAdoptions = [M.consumerAdoptions M2.consumerAdoptions];
	M.providerAdoptions  = [M.providerAdoptions  M2.providerAdoptions];
	M.providerExits      = [M.providerExits      M2.providerExits];
	M.laggedConsumers   = [M.laggedConsumers   M2.laggedConsumers];
	M.laggedProviders    = [M.laggedProviders    M2.laggedProviders];
	
	clear pred_aggVals M2;
end

%%% Load results from each counterfactual
cnfactLabels = cell(1,NumCnfacts);

for cidx = 1:NumCnfacts
	load(sprintf('%s/cnfacts/%s/cnfact_vals.mat', config.runFolder, cnfactNames{cidx}), 'cnfact_aggVals', 'cnfactName', 'cnfactLabel');
	cnfactLabels{cidx} = cnfactLabel;
	
	M2.consumerAdoptions = array2table(mean(cnfact_aggVals.consumerAdoptions, 2)); M2.consumerAdoptions.Properties.VariableNames = {cnfactName};
	M2.providerAdoptions  = array2table(mean(cnfact_aggVals.providerAdoptions,  2)); M2.providerAdoptions.Properties.VariableNames  = {cnfactName};
	M2.providerExits      = array2table(mean(cnfact_aggVals.providerExits,      2)); M2.providerExits.Properties.VariableNames      = {cnfactName};
	M2.laggedConsumers   = array2table(mean(cnfact_aggVals.laggedConsumers,   2)); M2.laggedConsumers.Properties.VariableNames   = {cnfactName};
	M2.laggedProviders    = array2table(mean(cnfact_aggVals.laggedProviders,    2)); M2.laggedProviders.Properties.VariableNames    = {cnfactName};
	
	M.consumerAdoptions = [M.consumerAdoptions M2.consumerAdoptions];
	M.providerAdoptions  = [M.providerAdoptions  M2.providerAdoptions];
	M.providerExits      = [M.providerExits      M2.providerExits];
	M.laggedConsumers   = [M.laggedConsumers   M2.laggedConsumers];
	M.laggedProviders    = [M.laggedProviders    M2.laggedProviders];
	
	clear cnfact_aggVals cnfactName cnfactLabel M2;
end

if includeBase
	myLabels = [{'Base case'} cnfactLabels];
else
	myLabels = cnfactLabels;
end


%%% Make time series plots
mytitle = '';

%% Plots in B&W for plot_seeding and plot9_10
plotTimeSeries(M.laggedConsumers, sprintf('%s/timeseries_laggedConsumers.pdf', outputFolder), mytitle, '', ...
	myLabels, {}, {'black', 'black', 'black'}, {}, 'northwest', [0, 20000], '', {'-', '--', ':'}, 50, 20, 15);
plotTimeSeries(M.laggedProviders, sprintf('%s/timeseries_laggedProviders.pdf', outputFolder), mytitle, '', ...
	myLabels, {}, {'black', 'black', 'black'}, {}, 'northwest', [0, 10000], '', {'-', '--', ':'}, 50, 20, 15);

end